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Statistical distribution of energy levels for Dirac fermions confined in a quantum dot is studied 
numerically on the examples of triangular and hexagonal graphene flakes with random electrostatic 
potential landscape. When increasing the disorder strength, level distribution evolves from Pois- 
sonian to Wigner, indicating the transition to quantum chaos. The unitary ensemble (with the 
twofold valley degeneracy) is observed for triangular flakes with zigzag or Klein edges and potential 
varying smoothly on the scale of atomic separation. For small number of edge defects, the unitary- 
to-orthogonal symmetry transition is found at zero magnetic field. For remaining systems, the 
orthogonal ensemble appears. These findings are rationalized by means of additive random-matrix 
models for the cases of weak and strong intervalley scattering of charge carriers in graphene. The 
influence of weak magnetic fields, as well as the strong-disorder-induced wavefunction localization, 
on the level distribution is also briefly discussed. 
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I. INTRODUCTION 

The advent of graphene pQ, an experimental model 
for two-dimensional massless Dirac fermions [2] , has pro- 
vided the unique opportunity the test the theoretical pre- 
dictions for such exotic particles from one perspective [3] , 
and to reexamine the classic effects of nanoscopic physics 
[4] from the other perspective. Theoretical predictions al- 
ready verified experimentally include Klein tunneling [5] 
pseudodiffusive shot noise U\ ; the universal quantum 
values of the conductivity [BHS] and the visible light opac- 
ity [9]. On the other hand, nanostructures in graphene 
show weak localization [10j . universal conductance fluc- 
tuations [TJJ , Aharonov-Bohm [T2TfT¥] and Josephson [To] 
effects, just to mention a few. 

These two perspectives unify in the issue of quantum 
chaos in Dirac billiards [16] , modelled experimentally 
within the graphene quantum-dot devices |17j . Early 
theoretical considerations [16] suggested that in Dirac 
billiards time-reversal symmetry (TRS) may be broken 
even in the absence of magnetic fields, leading to the 
unitary symmetry class. Existing Coulomb-blockade ex- 
periments [17] provide strong indications for quantum 
chaos in graphene, but without giving a clear identi- 
fication of the symmetry class. Several computational 
experiments were performed [18H20] showing, that level- 
spacing distributions for irregular or disordered graphene 
nanoflakes exhibit orthogonal symmetry class, as scatter- 
ing the carriers between K and K' valleys restores TRS 
[18j . In turn, when searching for the unitary symmetry, 
one should focus rather on open than closed nanosystems 
in graphene. 

In this paper we follow the line of approach established 
with Refs. [18H20] . but focus on highly-symmetric (tri- 
angular and hexagonal) graphene nanoflakes, in which 
transition to quantum chaos is driven by weak potential 
disorder [2TH24] attributed to the influence of substrate 
impurities (or ions). Our numerical results show, that 
albeit the orthogonal symmetry appears generically in 



closed graphene nanosystems, in a peculiar case of tri- 
anguiar nanoflakes with zigzag (or Klein 25j) edges and 
smooth impurity potential the unitary symmetry class is 
the relevant one. Also, in such a case the (approximate) 
twofold valley degeneracy is observed, as the scattering 
of carriers between the valleys is negligibly weak. When 
turning on the magnetic field, such a system transforms 
into a pair of two independent chaotic systems (one at 
each valley) each of which showing the unitary symmetry. 
On the other hand, edge defects at zero field increase the 
intervalley scattering, such that the twofold degeneracy is 
lifted up and TRS is restored (leading to the orthogonal 
symmetry class). These findings complement the dia- 
gram of possible transitions between symmetry classes of 
graphene nanoflakes (see Fig.[T]). 

The structure of the paper is as follows. In Section |H} 
we discuss the relevant symmetries of the Hamiltonian 
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FIG. 1: Transitions between symmetry classes and random 
matrix ensembles relevant for closed nanosystems in graphene 
characterized by the disorder strength, the intervalley scatter- 
ing rate, and (optionally) placed in the weak magnetic field 
B. Solid arrows in the right part indicate transitions already 
reported in the literature; dashed arrows in the central part 
indicate remaining transitions. 
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for Dirac fermions in graphene and present the tight- 
binding model for a potential disorder. This model is 
utilized in Section [III| to demonstrate the level repulsion 
appearing in graphene nanoflakes when increasing the 
disorder strength. Next, in Section |IV| we overview the 
basic random-matrix models for dynamical systems and 
apply them to describe the transition to quantum chaos 
in graphene nanoflakes of different shapes and bound- 
aries. In Section [VJ we investigate the two distinct tran- 
sitions between the orthogonal and the unitary symmetry 
appearing in triangular nanoflakes at zero or finite mag- 
netic fields, and compare the spectral statistics for each 
case with those obtained from relevant random-matrix 
models. The influence of disorder-induced localization 
on spectral statistics is discussed in Section [Vl| The con- 



of the Dirac equation "HDiracV" = EJtjj via 



elusions are given in Section VII 



II. 



DIRAC FERMIONS IN WEAKLY 
DISORDERED GRAPHENE 



The microscopic model of disorder in graphene 
nanoflakes, representing the random electrostatic po- 
tential landscape [HI [2lTf2"4"] . is presented in this Sec- 
tion. Depending on the disorder correlation length £, the 
model may represent the potential abruptly (£ -C a) or 
smoothly (£ 3> a) varying on the length-scale of the lat- 
tice spacing in graphene a — 0.246 nm. But first, let us 
briefly recall (after Ref. [13]) the discussion of possible 
symmetry classes of such nanosystems. 



A. Symmetries of the Hamiltonian 

The effective Hamiltonian for low-energy excitations 
of electrons in graphene in low magnetic fields |26j has 
a form of the Dirac Hamiltonian 



^Dirac 

VF(Px + eA x ) a x ®T z + v F (p y + eAy) a y ® t 

+ M(x,y)a z <g)T + U(x,y)a ® r , (I) 

where vp — 10 6 m/s is the energy-independent Fermi ve- 
locity, Pi — —ihdi (with i = 1,2) are in-plane momen- 
tum operator components, Oj and Tj (j = 1,2,3) are the 
Pauli matrices acting on sublattice and valley degrees 
of freedom (respectively), and cto (to) denotes the unit 
matrix. The electron charge is — e, the vector potential 
A = (A x ,A y ) defines perpendicular magnetic field via 
B z — e z ■ rot A = d x A y — d y A. x , whereas M(x,y) and 
U(x, y) are the mass term and the electrostatic poten- 
tial energy (respectively). The Hamiltonian ([I]) acts on 
spinors ij) = [ipA, ipB, V'A' V'b] T > wnere A/B is the the 
sublattice index and the primed and unprimed entries 
correspond to the two valleys. Two-component wave- 
function for a charge carrier in the position representa- 
tion *&(x, y) = [^a,^b] t is determined by the solution 



*(r) 



4>B 



iK'- 



(2) 



where r = (x,y) and K (K') stands for the position of 
K (K') valley in the momentum space. (We choose K = 
— K' = |f for the remaining parts of the paper.) Eq. 
([2]) allows one to discuss the position dependence of each 
spinor component of slowly varying on the scale of 
atomic separation a (apart from the full wavefunction 
^(r) varies abruptly on the scale of a). 

Symmetries of the Hamiltonian ([!]) are defined by the 
following antiunitary operations: standard time reversal 
T, and two "special time reversals" 

T ={o a ® Tx )C, (3) 
T s \ = -i(o- y ® t )C, T v = -i(o- <8> t v )C, (4) 

where C denotes complex conjugation. The mass term 
M(x, y)o- z ®TQ breaks the symplectic symmetry associated 
with Tsi, leading to the two distinct possible scenarios: 

(i) In the case of weak intervalley scattering, T v com- 
mutes with 'HDirac, so the system consists of two in- 
dependent subsystems (one for each valley). Each 
subsystem lacks TRS (even at zero magnetic field) , 
as T commutes only with full "HDirac- Because the 
Kramer's degeneracy (T^ = —I) [37], the Hamil- 
tonian of a chaotic system consists of two degen- 
erate blocks (one per each valley), each of which 
may be modelled by a random matrix belonging to 
the Gaussian Unitary Ensemble (GUE). The anal- 
ogous scenario was first considered by Berry and 
Mondragon [TB] for neutrino billiards, lacking the 
valley degree of freedom. When magnetic field is 
applied to the system, ifoirac no longer commutes 
with Tv and the valley-blocks are not degenerate. 

(ii) In the case of strong intervalley scattering caused by 
irregular and abrupt system edges (or by a poten- 
tial abruptly varying on the scale of atomic separa- 
tion) the two sublattices are nonequivalent, so both 
special time-reversal symmetries T s \ and T v became 
irrelevant. For B = \B Z \ — 0, T commutes with 
^Dirac leading to the orthogonal symmetry class 
and statistical properties following from the Gaus- 
sian Orthogonal Ensemble (GOE) of random ma- 
trices. When increasing B, transition GOE-GUE 
similar to that discussed earlier for Schrodinger sys- 
tems [5S] appears. 

The existing numerical studies for closed systems of ir- 
regular shapes [I8H20] show that the typical intervalley 
scattering time is always shorter than the time required 
to resolve a level spacing (Heisenberg's time) leading to 
the scenario (ii). The corresponding transitions between 
ensembles of random matrices are depicted in Fig. [l] with 
solid lines. Some features of the scenario (i) were found in 
open systems [151 HSj; for which the intervalley scattering 
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time needs to be compared with much shorter time char- 
acterizing the conductance (escape time). Such systems 
are, however, beyond the scope of this paper. We focus 
here on a peculiar case of regular and weakly-disordered 
nanosystems, for which the intervalley scattering itself 
may be strongly suppressed, providing the appropriate 
boundary conditions are chosen (see Appendix [A[) . 



B. Potential disorder in the tight-binding model on 
a honeycomb lattice 

The lattice Hamiltonian for disordered graphene in 
weak magnetic field reads 



H T BA = ^MA)|*)01+h.C.] 



{ij) 



+ J2 [Mv(r t ) + E/gatefc) + U imp (ri) (5) 



The complex hopping-matrix elements are given by 



-texp 



J J A • dr (with the flux quantum 



$0 = h/e ~ 4.14 x lCT 15 T-m 2 ) if the orbitals |i) and 
\j) are nearest neighbors on the honeycomb lattice (with 
t = ^V3hvp/a ~ 3eV), otherwise Uj = 0. (The symbol 
J2(ij) denotes that each pair (ij) is counted only once.) 
This represents a minimal form for the interaction be- 
tween lattice fermions and the magnetic fields (Peierls 
construction |29| ) in a framework of the tight-binding ap- 
proximation (TBA). The mass term is modelled within a 
staggered potential on a honeycomb lattice My(r,) [30] . 
which is positive (negative) if r 4 belongs to sublattice A 
(B). Typically, we put |My(rj)| < t if is the outermost 
atom position at zigzag edge, or My(rj) = otherwise. 
Such a simple choice was shown to reproduce the "infinite 
mass" boundary condition correctly for various scatter- 
ing problems [3T]. The physical origin of a staggered 
potential is usually related to the magnetic moments at 
the zigzag edges [55] • (Alternatively, high electrochem- 
ical potential of terminal atoms can also be attributed, 
for instance, to the hydrogen-edge passivation, see Ref. 
[33].) We also consider the case of My(rj) = at all 
lattice sites for the comparison. 

Finally, the electrostatic potential term in "Htba con- 
tains a contribution L/ ga t c from gate electrodes (slowly 
varying with the site position r^) and a random contri- 
bution U- lmp from impurities. For small nanoflakes one 
can choose Z7 gat0 ~ Uq = const, whereas a realization 
of disorder potential is generated by randomly choosing 
Ni m p lattice sites R„ (n = 1, . . . , N imp ) out of JV tot , and 
by randomly choosing the amplitudes U n € (—5, 5). The 
potential is then smoothed over a distance £ by convolu- 
tion with a Gaussian, namely 



^imp(r) = ^2 U n e *V 



2e 



(6) 



The special case of £ <C a, A^ rop = N to t corresponds 
to the Anderson model on a honeycomb lattice, consid- 
ered in Ref. [TS] on spectral statistics of graphene and 
nanotube-like structures. Earlier, the model constituted 
by Eqs. ( |5|6[ ) with £ ^> a was shown to reproduce basic 
transport properties of disordered mesoscopic graphene 
samples [TTJ [3T] [53] . Apart from a very recent work of 
Ref. [23], it has not been considered in the discussion of 
spectral statistics of nanoflakes so far. 

We further define the Fourier transform of two-point 
correlation function 



A 



Ntot Ntat 



{N tot hv F ) l=1 j=1 



xexp[iq-(ri-rj)], (7) 

where the system area A = ^y/3N tot a 2 , and the averag- 
ing takes place over possible realizations of the disorder 
^ (so (J7im P (r)) = 0). For the length scales large com- 
pared to £, the dimensionless correlator 



V3 Nut* /f j 2 
9 N tnt \ t ' 



., if £ < a, 

;^(£/a) 2 , if£»a, 



(8) 



becomes a representative measure of the disorder 
strength. For q ^ 0, we obtain Kq = Kq if £ -c a, or 
if q = K exp(— q 2 £ 2 ) if £ » a. (The latter justifies call- 
ing £ the 'disorder correlation length', as proposed earlier 
in the paper.) The numerical value of the ratio Kq/K a 
at q = ±K = (± I s , 0) roughly approximates the inter- 
valley scattering rate [23], and is as small as 2 x 10~ 6 
for £ = v3o (the value used for computer simulations 
presented in the remaining parts of the paper). 

Apart from negligibly weak intervalley scattering dis- 
cussed above, U- lm p also contributes to the mass term in 
^Dirac Q and thus breaks the symplectic symmetry as- 
sociated with t s i (|4) independently from the fact, that a 
similar effect may be caused by the system boundaries 
[34] . Namely, the effective mass term for low-energy ex- 
citations can be approximated by 

+ ^VC/i mp (fy)-(r i -r J ), (9) 

where and Vj are in-plane positions of atoms in the 
same unit cell [with i (j) belonging to the sublattice A 
(B)]. We further define r.y = (r^ + r 3 )/2. It is clear from 
Eq. |9|, that M e & ^ even for M v = 0. For higher 
energies, the symplectic symmetry is also broken by a 
nonlinear term appearing in the effective Hamiltonian de- 
rived from a nearest-neighbor tight-binding model [35j . 
Therefore, the structure of "Htba ^ provides additional 
reasons, for which energy levels of graphene nanoflakes in 



4 




TABLE I: Types and geometric sizes of nanosystems which 
spectral statistics are discussed in Sees. Ill V] 44 . 




-1.0 



FIG. 2: Systems studied numerically in the paper, (a), (b) 
Triangular and hexagonal nanoflakes with armchair, zigzag 
and Klein edges, (c) Typical impurity potential landscape 
(7i mp (r) for a triangular flake with armchair edges and 2106 
carbon atoms. The triangle height is H — 39 a ~ 10 nm, the 
impurity concentration and the disorder correlation length are 
N-aap/Ntot = 0.01 and £ = y/3a, respectively. 



the limit of quantum chaos, obtained numerically in the 
remaining parts of the paper, follow GOE or GUE statis- 
tics, depending whether charge carriers are scattered be- 
tween the valleys or not. 



C. Scope of the paper: Nanosystems considered 
and the numerical approach 

Graphene nanoflakes studied in the paper are shown 
schematically in Fig. [2] In general, we limit the discus- 
sion to the triangular and hexagonal flakes bounded en- 
tirely with armchair, zigzag, or Klein edges [see Figs. 
[2]Ja) and[2]Jb)]. Such a choice is related to the fact, that 
these three types of edges were observed using different 
microscopy techniques [55H55], In particular, STM mea- 
surements for quantum dots with well-defined edges and 
almost hexagonal shapes have been reported [39 . For 
this reasons, and because of numerous earlier theoretical 
works focused on graphene nanosystems with irregular 
edges [T8lI20| . it is worth to find out whether quantum 
chaos may even appear for highly-symmetric systems in 
the presence weak bulk disorder, and (if so) how symme- 
try classes of such systems are related to the boundary 
conditions? Later (in Sec. [V]) we extend our analysis on 
flakes with some randomly-distributed edge vacancies, to 
make link with the results of Refs. 



Shape 


Edges 


Height 


# of atoms 


Ave. size vA 






H/a 


ATtot 


xa 1 


xnm 1 


Triangle 


armchair 


78 


8268 


59.8 


14.7 






156 


32760 


119.1 


29.3 




zigzag 


45^3 


8278 


59.9 


14.7 






90^3 


32758 


119.1 


29.3 




Klein 


45.5v/3 


8548 


60.8 


15.0 






90.5^3 


33298 


120.1 


29.5 


Hexagon 


armchair 


64.5 


8322 


60.0 


14.8 






130.5 


34062 


121.4 


29.9 




zigzag 


42^3 


10584 


67.7 


16.7 






84^ 


42336 


135.4 


33.3 




Klein 


41.5^3 


10332 


66.9 


16.5 



As described in Sec. II B weak bulk disorder in our 



nanosystems is introduced within random electrostatic 
potential landscape. The potential fluctuations, usually 
attributed to the influence of substrate impurities, give 
origin to the so-called "puddles," i.e., spatial fluctuations 
in carrier density observed in numerous experiments |40F 
W2\ . We note here, that similar fluctuations may also 
arise from out-of-plane lattice deformations, modifying 
the electrostatic potential term [43] . 

An example of the potential given by Eq. ([6]) is shown 
in Fig. [2jc) . For demonstrating purposes, we took a rel- 
atively small triangular flake with armchair edges, which 
consists of iV tot = 2106 carbon atoms, corresponding to 
the triangle height H — 39 a ~ 10 nm. [When analyzing 
statistical distributions of energy levels, we choose the 
systems significantly larger, see Table [I]] The remaining 
parameters are the impurity concentration N- lmp /Nt t = 
0.01 and the disorder correlation length £ = \^3a, corre- 
sponding to Kq ~ 1.16 (5/t) 2 . Although the impurities 
visualised in Fig. [2][c) are relatively well-separated from 
each other, as well as K <C 1 for a typical value of 



5/t = 0.1 used in the simulations, we show in Sees. Ill 
and |IV| that such a weak disorder may lead to clear signa- 
tures of quantum chaos, providing the system considered 
is sufficiently large. 

Linear sizes of nanosystems listed in Table [I] are close 
to these reported in Ref. [T7J. Such systems contain 
10 4 < iVtot < 10 5 carbon atoms, making possible to re- 
produce, on a finite lattice, several features of a con- 
tinuous system (described by the Dirac theory) with 
a good accuracy [HJ OH]- On the other hand, the val- 
ues of ATtot 10 4 combined with the presence of a ran- 
dom potential landscape £/i mp (r) make it difficult 
(unless impossible) either to utilize ab-initio methods 
for carbon-based nanosystems [32 |45l |46] or to employ 
semiclassical theory for generic Dirac billiards in uni- 
form potentials presented in Ref. [37]. For these rea- 
sons, our method of approach is founded on a numeri- 
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cal diagonalization of tight-binding Hamiltonians Wtba 
([5]) for different nanosystems and different Ui mp (r). In 
brief, when analysing the spectral statistics of smaller 
systems (N to t J$ 10 4 ) we took 200 — 400 randomly-chosen 
CAmp(r) for each disorder strength quantified by the cor- 
relator Kq. For larger systems (iVtot S> 10 4 ) we took 
just one J7i mp (r) for each K , essentially reproducing the 
experimental situation of Ref . |17) , where the Coulomb- 
blockade spectrum of a single device was obtained. This 
allows us to verify, whether spectral statistics obtained 
for the ensemble of smaller systems coincide with spec- 
tral statistics obtained for a single, however much larger, 
system of each kind. 

More details on our numerical approach and methods 
of data analysis are provided in Sec. |IV| But first, let 
us briefly discuss a level structure of integrable Dirac 
systems on the example of perfect triangular graphene 
nanoflakes, and demonstrate how the level-repulsion ap- 
pears when weak disorder is included. 



III. LEVEL REPULSION IN TRIANGULAR 
GRAPHENE NANOFLAKES 



Energy levels of perfect triangular nanoflakes 
and the effect of weak disorder 



Energy levels of triangular nanoflake with armchair 
edges were recently found analytically, in the absence of 
disorder, by Rozhkov and Nori [48]. Close to the Dirac 
point, exact energies can by approximated as [461 148j 





Poisson 




< 





;gue/goe; 



rp± i^Dirac 

m,n — m,n,± 



±- 



2nt 



vOTot 

with 1 ^ m ^ 11. 



mn, (10) 
(11) 



-E° n jfj_ corresponds to eigenenergies of Dirac particles 
in a triangular cavity with upper (lower) sign valid for 
electrons (holes). For zigzag edges, an analytic solution 
of a tight-binding model is missing, but an approximation 
( 10 1 remains valid with (m, n) obeying the conditions 



777 ^ 1 and 77 ^ 2777. 



(12) 



The above also applies for Klein edges. 

Numerical examples of energy levels for triangular 
Dirac cavities and graphene nanoflakes are provided and 
discussed in Appendix [B] Here we only mention, that for 
a large system vast majority of such energy levels shows 
the fourfold degeneracy (per each direction of spin) in the 
absence of disorder and for My = 0. For zigzag (or Klein) 
edges, this degeneracy can be easily attributed to special 
time-reversal symmetries l~ s \ and T v Q. For armchair 
edges, the symmetry associated with T v no longer ap- 
plies, but the new twofold degeneracy E^™ c — E^"^ n 
appears instead. This is an accidental degeneracy, which 
is lifted if one includes the disorder. Also, the degeneracy 
associated with 7~ s \ is lifted in the presence of disorder for 
all types of edges, as the disorder potential leads to the 



0.3 




0.2 < 



FIG. 3: Evolution of energy levels for a triangular flake with 
armchair (left) and zigzag (right) edges containing iVtot = 
8268 and 8278 carbon atoms (respectively) when varying the 
disorder amplitude 8 (see Ref. A9 ). The index values ( m, n ) 
refers to electronic energies for Dirac cavities £ l m,»v4- dlOp- 
The flow-chart diagram (top) illustrates the transitions to cor- 
responding chaotic ensembles (see Fig. [TJ. 



effective mass term ([9|. In turn, only the twofold val- 
ley degeneracy (associated with %) for zigzag or Klein 
edges appears to be robust against weak potential disor- 
der. Such a degeneracy is also insesitive to the staggered 
potential My 7^ appearing at zigzag edges. 

The behavior of energy levels for triangular flakes with 
gradually increasing disorder is illustrated in Fig. [3j We 
took triangles containing N to t = 8268 and 8278 atoms 
(with armchair and zigzag edges, respectively), a sin- 
gle disorder realization defined by positions of impuri- 
ties R„ and their relative amplitudes U n /S E (—1, 1) for 
1 ^ 77 ^ -/V; mp ([6]), and varied the absolute amplitude 5 
|49) . Such a procedure may correspond to varying the 
distance between a graphene sample and its substrate in 
the case of free-standing samples [50], or to modifying 
the charge screening efectiveness in recently fabricated 
graphene- hBC heterostructures [51] . 

It is clear from Fig. [3] that due to low density of states 
in the vicinity of the Dirac point, Eq. ( 10 ) describes only 
the few lowest-lying states accurately for iVtot < 10 4 (no- 



tice the index values (m, 77) provided in the central part 
of a plot). The basic features of the level structure in 
the presence of disorder are, however, correctly repro- 
duced. For armchair edges, degeneracies discussed above 
are lifted for a relatively weak disorder strength, and one 
can directly compare a level sequence obtained from the 
numerical diagonalization of Htba (|5J) with predictions 
of the random matrix theory (RMTjiFor zigzag edges, 
the twofold (approximate) valley degeneracy of each level 
remains even for stronger disorder (unless \E\ < 8, see 
dashed line in the right panel of Fig. and the additional 
effort in the data analysis is required. Also, when increas- 
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ing the disorder strength, we observe avoided crossings 
signaling the level repulsion, characteristic for chaotic 
quantum systems [57]- Moreover, the repulsion is notice- 
ably stronger for zigzag edges and |J5| > S than for the 
other cases. The corresponding sequence of energy levels 
for Klein edges is not shown in Fig. [3j but it evolves with 
the increasing disorder strength in an identical manner 
as the level sequence for zigzag edges. 

These observations, together with the discussion of 
symmetry classes provided in Sec. [TTJ suggest that differ- 
ent ensembles of random matrices are capable of describ- 
ing level sequences such as depicted in Fig. [3j Namely, 
GUE with the approximate valley degeneracy shall ap- 
ply for triangular flakes with zigzag (or Klein) edges in 
the energy range \E\ > S, whereas GOE (without degen- 
eracy) shall apply otherwise. Resulting transitions be- 
tween the Poissonian and different Gaussian ensembles 
expected when increasing the disorder strength are pre- 
sented in a flow-chart diagram at the top of Fig. [3] We 
further verify these predictions in Sec. IV with the help 
of statistical analysis of the level sequence. 



B. Level-spacing distribution 

Basic statistic which distinguishes the spectra of inte- 
grable systems from chaotic ones is the level-spacing dis- 
tribution P( k ^(S). By definition, p( k \S)dS represents 
the probability, that the quantity (p(E)) (Ei+k — Ej) / k 
is located in the interval (S, S + dS), where Ei+k—Ei 
is the distance between fc-th neighbors in the level se- 
quence Ei ^ Ei ^ . . . , and (p(E)) is the average density 
of levels in the energy interval (E, E + dE) . Unlike for 
Schrodinger systems, such as the two-dimensional elec- 
tron gas (2DEG) where (p(E)) is usually assumed to be 
constant, in bulk graphene 

(p(E)) ~ p hulk (E) = -—^—\E\. (13) 

7T (hVF) 

Due to a quasirandom character of the Hamiltonian, sta- 
tistical properties of energy levels of large quantum sys- 
tems are described by RMT [53]. For instance, generic 
integrable systems are described by the Poisson distribu- 
tion, namely [54] 



P {k) (S) 



5 fc " 1 exp(-fcS*), 



(fc-1)! 



(14) 



with k = 1,2,.... For instance, the nearest-neighbor 
spacing distribution P^(S) = Pp i(S) — exp(— S). In 
contrast, nearest-neighbor level spacings P^(S) of clas- 
sically chaotic systems which preserve (or break) TRS 
may be approximated by the so-called Wigner surmise 
for GOE (or GUE) of random matrices 



Pgoe(S) = -Sexp 



GUE 



(S) 



32 

^2 



S^exp 



irS 2 
4 

AS 2 



(15) 
(16) 




0.8 




12 3 

- — I 1 1 h-^—i 1 h 

(c) TBA , ... ... "r (f) TBA .. 
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Poisson 
GUE 
, t = 32758 
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FIG. 4: Integrated level-spacing distributions C"' 2 '(S) for 
triangular flakes with armchair (a)-(c) and zigzag (d)-(f) 
edges, both having th e ar ea of A — (120 a) 2 [see Table I . 
(a),(d) Levels E^n% (10' f° r Dirac cavities. (b),(e) Levels 
obtained from diagonalization of Htba ^ in the absence of 
disorder (Kq = Q); and (c),(f) for infinitesimally-weak disor- 
der (K ~ 10~ 3 ) |52| . Insets show level-spacing distributions 
P( 1,2 \S). Numerical results are shown with solid black lines. 
Remaining lines are for Poisson distribution (blue dashed) 
and for the Wigner surmise for GOE or GUE (blue dotted). 



One notes that Pp i(S) — 1 — S for S — > and is max- 
imal for S = 0, i.e., integrable systems exhibit level at- 
traction. In contrast, Pqoe{S) oc S for S — > showing 
linear level repulsion, whereas Pgue(S) oc S 2 for S — >• 
showing stronger (quadratic) level repulsion than GOE. 
The above holds true for a simple sequence of energy lev- 
els [53], i.e., a sequence in which all levels have the same 
values of quantum numbers corresponding to strictly con- 
served quantities (resulting from the symmetry of the 
system). As discussed earlier in this Section, the valley 
index in graphene needs to be regarded as an example 
of such a quantity for triangular nanoflakes with zigzag 
or Klein edges even in the presence of disorder (see also 
Appendix [b]) . 

Before analyzing the level-spacing distribution as 
a function of disorder, we first briefly discuss that charac- 
terizing perfect (or almost perfect) triangular nanoflakes 
in graphene. In principle, equilateral triangles are not 
examples of generic integrable systems as they show 
number-theoretic degeneracies |55| . As shown in Ap- 
pendix [C] the average degeneracy of levels given by Eq. 
( 10 1 in a finite energy interval is divergent for A^ot ~> oo 
(a phenomena known as level clustering) and P^ (S) 
does not exists (as for similar corresponding Schrodinger 
systems [515]). For a large but finite N to t, the values of S 
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occurring for E^ 1 ™ ± are equal to [57] 



Sq ~ 18 



q~ 0.302 q, with q = 0,1,2, 



As a result, the integrated spacing distributions 



C( fe )(S*) = / P (k) (S')dS' 



(17) 



(18) 



show abrupt steps at S = S q . 

Such steps are clearly visible in Figs. |4^a) andQd), 
where we plot C W (5) with k = 1 (fc = 2) for a finite 
Dirac cavity with armchair (zigzag) boundary conditions 
[solid lines]. We took all energy levels given by Eq. (10) 
with \E%™ ± \ < t/2 for Wot = 32760 (or Wot = 32758) 
and used ptmik ( 13 ) to unfold the spectrum. The statistic 
CW(S) is used in case of zigzag boundary conditions due 
to the twofold valley degeneracy of each electronic level. 
The steps of (S) are followed by peaks of p(*0 (S) 
at S = S q (see insets). The discrete structure of spac- 
ing distributions gets smeared out when considering elec- 
tronic levels of %tba ([5| even in the absence of disorder 
(see Figs.Qb) andQe)). This is because %tba leads to 
the nonlinear dispersion relation and number-theoretic 
degeneracies no longer apply for triangular nanoflakes 
(see Appendix[C| . For this reason, such highly-symmetric 
nanosystems in graphene appear to be generic integrable 
systems, with Poissonian distribution of level spacings, 
providing the degeneracies associated with special time- 
reversal symmetries Q are properly taken into account. 
In fact, only small deviations from Pp i{S) [dashed lines] 
are visible due to finite numbers of energy levels consid- 
ered (see insets in Figs.Wb) and[4je)). Finally, in Figs. 
4[c) and gf) we plot C^(S) for triangular nanoflakes 
with an innnitesimally-weak bulk disorder (Kq ~ 10 ) 
[52"] to illustrate the level repulsion for S -C 1. The re- 
pulsion is noticeably stronger for a triangle with zigzag 
edges than for a triangle with armchair edges, suggesting 
that different symmetry classes apply in these two cases 
(see spacing distributions obtained from the Wigner sur- 
mise for GOE and GUE; dotted lines). The evolution of 
spacing distributions with increasing disorder strength is 
analyzed in a quantitative manner in Sec. |IV[ 



IV. RANDOM MATRICES AND SPECTRAL 
STATISTICS OF DISORDERED SYSTEMS 

This Section and Sec. W\ present the central results of 
the paper. We start from a brief description of basic 
additive random-matrix models [27J [55] capable of repro- 
ducing the evolution of spectral statistics when dynamic 
system undergoes transition to quantum chaos or transi- 
tion between different symmetry classes. Next, we apply 
these models to parametrize the transition to quantum 
chaos in weakly-disordered graphene nanoflakes. 



A. Additive random-matrix models and transitions 
between ensembles 

When generic integrable system undergoes the tran- 
sition to quantum chaos, its spectral properties may be 
reproduced by the following random Hamiltonian 



H{\) = 



H° + XV 



(19) 



where H is diagonal random matrix, which elements 
follow a Gaussian distribution with zero mean and the 
variance ((i/?) 2 ) = 5ij, the parameter A <E [0,oo], and 

V = is a member of one of the Gaussian ensembles. 
In particular, for the transition Poisson-GOE, elements of 

V are real numbers chosen to follow a Gaussian distribu- 
tion with zero mean and the variance (V?) = (1+Sij)/N, 
where N is the matrix size. Analogously, for the tran- 
sition Poisson-GUE, elements of V are complex num- 
bers which real and imaginary parts are generated inde- 
pendently according to Gaussian distribution with zero 
mean and the variance ((ReVy) 2 ) = (1 + Sij)/2N and 
((ImVy) 2 ) = (1 - S tj )/2N, respectively [51]. 

For N — 2, the nearest-neighbor spacing distribution 
for the Hamiltonian ( 19 ) can be found analytically and 
reads, for the transition Poisson-GOE I5JH, 



-Ppoi-GOE(A; S) 



~u{\) 2 S~ 




' u(X) 2 S 2 ' 


A 


exp 


4A 2 



< / drj exp(-?7 2 — 2Xr])I Q 
Jo 



rju(X)S 
A 



(20) 



Iq(x) is the modified Bessel function of the first kind; 
u(X) = tJttU{— 1, 0, A 2 ) with U(a,b,x) the confluent hy- 
pergeometric function [BO]- For the transition Poisson- 
GUE we have [58] 



~a(X) 2 S~ 




' a(A) 2 S 2 " 


A 


exp 


2A 2 



-Ppoi-gue(A; S) — 



3 exp I — Xij — 

x / dr/ — sinh 

/o V 



rja(\)S 
A 



(21) 



where the coefficient a(A) is expressed by the error func- 
tion ed(x) = (2/y/W) J*exp(-t 2 )dt as 



a(X) 



exp 



1-erf 



V2 



+ A 2 



' dx exp(-Ax) exf 
X 



(22) 



In particular, for A = Eqs. (20) and (21) both restore 
the Poissonian distribution Pp j(5) = exp(— S). For 
the opposite limit (A — > oo) we have Ppoi-GOE(S') 
Pgoe(S) 



|15|) and Ppo^gue(S) ~ P G ue(S) (16), re- 
producing the Wigner surmise for GOE and GUE ma- 
trices (respectively). For < A < oo, Eq. p0| de- 
scribes level-spacing distributions interpolating between 



Poisson and GOE statistics, with Pp i-GOE(A; S) oc S/X 
if S < X < 1, or P{X-S) oc S if S < 1 < A. Analo- 
gously, Eq. ( plj ) describes level-spacing distributions in- 
terpolating between Poisson and GUE statistics, with 
Ppoi-GUE (A; S) oc S 2 /X if S < X < 1, or P(A; 5) cx S 2 
if S" <C 1 < A. In principle, for any A > the distribu- 
tions p0| and (21) both exhibit qualitatively the same 
level repulsion (i.e., linear or quadratic) as the Wigner 
surmise for corresponding Gaussian ensembles. 

For a sake of completeness, we also mention that the 
random-matrix model of the form given by Eq. (191 but 
describing the transition GOE-GUE (i.e., for H° a mem- 
ber of GOE, V a member of GUE, and N = 2) gives 
a simple expression for the nearest-neighbor spacing dis- 
tribution \M 



GOE-GUE 



(A; S) 



A 2 



Sc 2 (A) 



x exp 



S 2 c 2 {X) 



erf 



Sc{X) 



(23) 



with 



=(A) = 



tt(2+A 2 ) 



\ 2 ( ( A 

1 arctan — = 

t V VV2 



^2A 
2 + A 2 



(24) 



By varying the parameter A € (0, oo) one gets a family 
of distributions interpolating between Wigner surmises 
Pgoe(S) ^ and P GUE (S) pj). 

Despite Eqs. @, ((2TJ, and ^23|) are exact for 2x2 
random matrices only, they can also be utilized to 
parametrize transitions between ensembles of large ran- 
dom matrices. It was show numerically, that Px(X&t;S) 
with Afit oc VNX and X = Poi - GOE, Poi - GUE, or 
GOE — GUE, provides approximations of the nearest- 
neighbor spacing distributions of random matrices given 
by Eq. ( 19 ) for N 3> 1 with an astonishing accuracy [55] , 



Moreover, such approximations were applied to describe 
the energy spectra of various dynamic systems undergo- 
ing transitions between symmetry classes [371 HH liH • 
In the remaining part of the paper we show, that dis- 
tributions Px(X&t',S) are also relevant when discussing 
transitions between symmetry classes for Dirac fermions 
confined in graphene nanoflakcs. 



B. Energy-level distributions and transition to 
quantum chaos in graphene nanoflakes 

The evolution of level-spacin g distributions pW(5) ( or 
pW{S)) with the increasing disorder strength (quantified 
by K ) is illustrated in Fig. [5] [64] on the example of 
a triangular nanoflake with armchair (or zigzag) edges 
containing N tot — 32760 (or N tot = 32758) atoms (see 
Table |l| . To unfold the spectra of finite systems in the 
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FIG. 5: Level-spacing distributions P^ 1,2 '(S) for triangu- 
lar nanoflakes with armchair (a)-(c) and zigzag (d)-(f) edges. 
The flake area is A — (120 a) 2 . The disorder strength Kq is 
varied between the panels [64 . Numerical results are shown 
with black solid lines. Red solid lines show t he b est-fitted ap- 
proximating distributions -Ppoi-GOE(Afi t ; S) (20 1 [panels (a)- 
(c)] or Ppoi-GUE(A fi t;S') [panels (d)-(f); 



fied for each plot. The ot 



with Afit speci- 
ler lines are same as in Fig. [4] 



presence of disorder, we use an approximating formula 
for the average density of states 



(p(E)) =s Po 



A. 



■off 

~A 



Pbulk(P), 



(25) 



with pbuik(P) given by Eq. (13 1. The constant term po 
and the effective flake area A c s !$ A are determined via 



least-square fitting of Eq. (25) to the actual (p(E)) ob- 
tained numerically for a particular realization of U lmp (r) . 



For the energy range considered, Eq. ( 25 ) provides a rea- 



sonable approximation of (p(E)) obtained for disordered 
graphene with various approaches, including analytical 
calculations employing Born approximation [65) or STM 
measurements for epitaxial graphene samples |66) . 

The presentation in Fig. [5] starts for infinitesimally- 
weak disorder K ~ 10~ 3 [panels (a) and (d)], same for 
which integrated spacing distributions C' 1,2 '(5} are plot 
in Figs, gc) andfflf). Similarly as for C^' 2 >(S), the 
values of p( 1,2 )(S) [black solid lines in Fig. [H] are inter- 
mediate between these corresponding to the Poisson dis- 
tribution [blue dashed lines] and to the Wigner surmise 
for GOE or GUE [blue dotted lines]. The level repulsion 
is clearly visible at S <C 1 for both systems, but signifi- 
cantly stronger for the triangle with zigzag edges than for 
the triangle with armchair edges. Next, we enlarge the 
disorder strength Ko by factor 4 between the consecutive 
panels in Fig. [5j (a)-(c) for armchair edges, and (d)-(f) 
for zigzag edges. The level-spacing distributions converge 
to the corresponding Wigner surmises for Kq > 0.01, 
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FIG. 6: Least-squares fitted parameters Afi t for (a) transi- 
tion Poisson-GOE and (b) transition Poisson-GUE [see Eqs. 
(20 1 and (21 1, respectively] as functions of disorder strength 
for systems of Table [T] [67 . Datapoints in panel (a) corre- 
spond to triangles with armchair edges (A iVtot = 8268, 
A^tot = 32760) and hexagons with armchair (O Ntot = 8322, 
• iVtot = 34062), zigzag (□ iV tot = 10584, N tot = 42336) 
and Klein (0 A" to t = 10332) edges. Datapoints in panel 
(b) correspond to triangles with zigzag (A A*t ot = 8278, A 
Ntot = 32758) and Klein (V N tot = 8548, T N to t = 33298) 
edges. Logarithmic scales are used on all axes. Lines repre- 
sent the best fitted power-law relations (see Table |H| . 



TABLE II: Least-square fitted parameters in Eq. ([26j) cor- 
responding to lines in Figs. [6ja) and|6|b). The values for 
hexagons with Anderson-type disorder (£ = 0, N mp — N tot ) 
are provided in the last row for comparison. Numbers in 
parenthesis are standard deviations for the last digit. 



Shape 


Edges 


Xi 


a 


Triangle /hexagon 
Hexagon 
Triangle 


armchair 
zigzag/Klein 
zigzag/Klein 


0.033(4)° 
0.018(3) a 
0.082(8) b 


0.59(3) a 
0.60(3) a 
0.43(2) 6 


Hexagon, £ = 


zigzag 


0.046(3)° 


0.59(l) a 



"Transition Poisson-GOE. Transition Poisson-GUE. 



see Figs. IBTc) and [5^f ). The interpolating distributions 
-Ppoi-GOE(Afit;S') pOl and -Pp i-GUE(Afi t ; S) plj) with 



Afit obtained by least-square fitting [red solid lines] pro- 
vide good approximations of the actual P^' 2 '(S) for all 
values of Kq. In particular, the character of level re- 
pulsion for small S, which is approximately linear for 
transition Poisson-GOE and approximately quadratic for 
transition Poisson-GUE, is well-reproduced with the nu- 
merical data for triangles with armchair and zigzag edges 
(respectively). These are the numerical evidences show- 
ing, that symmetry classes of weakly-disordered trian- 
gular nanoflakes in graphene remain the same as pre- 
dicted for chaotic Dirac billiards with appropriate bound- 
ary conditions (see Sec. [Il] and Appendix|A|, namely: the 
orthogonal symmetry applies for armchair edges or the 
unitary symmetry applies for zigzag edges. (For the lat- 
ter case, the symmetry class is also insensitive to the 
staggered potential My ^ at the system boundary.) 
Similar agreement between -Pp i-GOE(Afit; S) (or 



-Ppoi-GUE(Afit; S)) and the actual level-spacing distribu- 
tions was observed for all nanosystems listed in Table [I] 
In Fig. [6] we plot the values of Afi t as functions of the total 
disorder strength N tot Ko [57]. We find such an extensive 
quantity makes it possible to identify the approximating 
power-law relations 



Afit — Afit(C) = AiC" with ( = N tot K . 



(26) 



The coefficients Ai and a (provided in Table |n|) still dif- 
fer between the systems with different shapes or bound- 
ary conditions, but remain unchanged when varying N tot 
and Kq independently with the remaining parameters 
fixed. Surprisingly, the datapoints obtained for the en- 
tire collection of quantum billiards listed in Table|l]group 
along just three distinct lines on the log-log plot (see 
Fig. [6]). The datapoints corresponding to either triangles 
or hexagons with armchair edges (showing the transi- 
tion Poisson-GOE when increasing the disorder strength) 
may be approximated by Afi t (C) (26) with the param- 
eters Ai and a given in the first row of Table [TT] [red 
solid line in Fig. [6]ja)]. The datapoints corresponding 
to hexagons with zigzag or Klein edges (also showing the 
transition Poisson-GOE) may be approximated by Agt(C) 
p6| with Ai and a given in the second row of Table |B] 
[red dashed line in Fig. [6f a)] . We notice, that the best- 
fitted values of the exponent a for these two situations 
equal to apoi-GOE — 0.6 within the range of errorbars. 
Also, the same value of a was obtained for hexagons with 
Anderson-type disorder (see the last row in Table |ll|), sug- 
gesting that it is specific for nanosystems undergoing the 
transition Poisson-GOE, regardless the microscopic de- 
tails of the disorder realization. Finally, for triangles with 
zigzag or Klein edges (showing the transition Poisson- 
GUE) the corresponding datapoints shown in Fig. [6^b) 
may be approximated by Afit(C ) PQf with parameters 
given in the third row of Table |II| [red solid line] . The 
exponent a is equal to ap i_GUE — 0.4 in this case. 

The numerical results presented in Figs. [5] and [6] con- 
stitute an onset of transition to quantum chaos in highly- 
symmetric graphene nanoflakes with a weak potential 
disorder. The actual level-spacing distributions follow 
-Ppoi-GOE(A; S) 



(20) 



or -Ppoi-GUE(A; 5) (21) obtained 
from basic random-matrix models, which are applicable 
for generic quantum system in the orthogonal or the uni- 
tary symmetry class. Depending whether the interval- 
ley scattering is strong or weak in a particular graphene 
nanosystem, we have P^(S) ~ fp i-GOE(Afit; S) or 
P( 2 '(S) ~ -Ppoi-GUE(Afit; S) (the systems conserving val- 
ley pseudospin show an approximate twofold degener- 
acy of each energy level even in the presence of disor- 
der). The parameter Afit oc (N tot K ) a , with the expo- 
nent a taking one of the two values: apoi-GOE — 0.6, or 
apoi-GUE — 0.4. These characteristics of transition to 
quantum chaos arc further supported by the behavior of 
more distant spacings distributions briefly discussed in 
the next subsection. 
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C. Spectral rigidity 

A customary measure of spectral fluctuations on scales 
longer than spacings described by P^(S) or p( 2 \S) is 
provided by the spectral rigidity A 3 (L) defined by Dyson 
and Mehta ESI 



A 3 (L) = \ ( Min [ dx [Af(x +x)-ax-bf 

L \(a,b) J-L/2 i 



(27) 



where x = (AT(E)) and Af(E) denotes the number of 
energy levels having energy between E m - ln > and 
B max > E > -Emin- In turn, the average (Af(E)) = 
dE(p(E)} with (p(E)} approximated by Eq. (25 1. 

[For negative E, (N(E)) = - f E Emin dE(p(E)) .} The 
spectral rigidity A 3 (L) represents the- least square devia- 
tion of the actual spectral staircase N(E) from the best- 
fitting function ax+b over a range x € (xq—L/2, xq+L/2) 
(the averaging in Eq. (27 1 runs over the interval center 
x ). Theoretical expectations for A 3 (L) are [681169] 



A {X) (L) 



L 
15 



— l —r [ dS(L- Sf(2L 2 -9LS -3S 2 )Y X (S), with X — Poi, GOE, GUE; 
15i 4 Jo 



Y Poi (S)=0, Y GO e(S) = 



sin(7r>S') 


2 d 


sin (ttS) 




irS 


+ dS 


ttS 


f 

Js 



sin (nt) 



nt 



dt, y GUE (s) = 



sin(vrS') 
ttS 



(28) 
(29) 
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FIG. 7: Spectral rigidity Aa(L) for triangular nanoflakes 
with armchair (a) and zigzag (b) edges. Datapoints show 
the results obtained numerically using the same system pa- 
rameters as in Fig. [5] Lines correspond to the theoreti- 
cal expectations [see Eqs. ( |28|29[ )] for Poisson distribution 
(black solid), GOE (red dashed), and GUE (red dotted). 
(Notice that theoretical curves on panel (b) are rescaled via 
A^ X \L) = 4A^ x \L/2) due to the twofold valley degeneracy 
of each energy level.) 



In particular, for the Poisson distribution we have 



(L) 



For X = GOE or GUE, A { 3 X) (L) 



Y^L for L <C 1 [7U], and the exact values for larger L can 



be obtained numerically from Eqs. ( 28|29 |. 

Figs. [7]Ja) and [7jb) show the spectral rigidity A 3 (i) 
for graphene nanoflakes same as in Figs. [5Va) and [5jb) 
(see also Ref. [53]). It is clear from Fig. HTa), that the 
spectral rigidity for triangles with armchair edges gradu- 
ally evolves, with increasing K 0l from the straight line for 
Poisson distribution (black solid line) to the curve depict- 
ing the theoretical prediction for GOE (red dashed line) 
obtained from Eqs. ( |28|29 ). Some deviations visible for 
L > 10 can be attributed to a finite system size (see also 
the second paper in Ref. [5U] ). For triangles with zigzag 



edges, the approximate twofold valley degeneracy is ob- 
served, and the theoretical predictions A 3 (L) need to 

be replaced by A 3 (L) = 4A { 3 X) (L/2), drawn in Fig.^b) 
for X = Poi (black solid line) and X = GUE (red dot- 
ted line). The evolution of the actual A 3 (L) between 
the limiting theoretical curves is observed also in these 
case, showing the convergence to the predictions for GUE 
is reached, in the range L < 10, for as small disorder 
strengths as K ~ 10~ 2 . 



V. TRANSITION GOE-GUE IN TRIANGULAR 
GRAPHENE NANOFLAKES 

In this Section, we first discuss the transition GOE- 
GUE at zero magnetic field on the example of a triangular 
graphene nanoflake with zigzag edges, a finite number of 
edge vacancies N vac (which modifies the intervalley scat- 
tering rate) , and the bulk disorder strong enough to drive 
the system into chaotic regime. Then, we demonstrate 
the evolution of spectral statistics with the increasing 
magnetic flux $ piercing the system. 



A. Level-spacing distributions in the presence of 
edge vacancies 

Spectral characteristics of triangular nanoflakes with 
zigzag edges and a finite number (iV vac ) of vacancies, 
randomly-distributed at the system boundary, are pre- 
sented in Figs. [8] and [9] In particular, the evolution 
of nearest-neighbor spacing distribution PW(S) for < 
Avac ^ 10 is shown in Figs. |8ja)-(c), where we have cho- 
sen the total disorder strength iVtot-fth — 1-3 x 10 3 [71]. 

Earlier, we found that second-neighbor spacing distri- 
bution p( 2 \S) for the system without vacancies (iV vac = 
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(d) N vac =0 ., 

i >. flt = 0.051 




FIG. 8: Level-spacing distributions P^\S) (a)-(c) and 
P^(S) (d)-(f) for triangular nanoflakes with zigzag edges, 
the area A = (120a) 2 , and the disorder strength fixed at 
K ~ 0.04 (corresponding to N tot K ~ 1.3 x 10 3 ) [71] • The 
number of edge vacancies iVvac is varied between the panels. 
Numerical results are shown with black solid lines. Th e ot her 



lines correspond to empirical distributions pl x I(A; S 1 ) (32 1 [or 



pJ?^(A; S) ( 33 1] with least-square fitted A = \a t (red solid), 
A = (blue aashed), and A = oo (blue dotted). 



0) can be approximated by Pp i-GUE(Afit, S) pi] ) with 
Afit = 1-62 (see Fig. [6jb) ) , truncating the distribution 
following from GUE of random matrices with approxi- 
mate twofold (valley) degeneracy of each level. This ob- 
servation is further supported by the bimodal structure 
of PW(S') visible in Fig. ^a). The first mode of the 
distribution obtained numerically using levels with ener- 
gies 0.1 ^ \E\/t ^ 0.5 (black solid line) is centered at 
S ~ 0, and corresponds to the contribution of odd spac- 
ings, separating two almost-degenerated copies of each 
energy level belonging to K and K' valleys. (Notice, that 
weak intervalley scattering is present in a finite lattice 
system even for iVvac = 0.) The second mode, centered 
at S ~ 2, corresponds to the contribution of even spac- 
ings, and reproduces the Wigner surmise for GUE (l6| ) 
scaled according to P^(S) ~ \P GVE (S/2) (blue dashed 
line) with an excellent accuracy for S > 1. For N vllc > 0, 
the contribution of odd spacings gets smeared out and 
the distribution P^(S) follows the corresponding GOE 
statistic (blue dotted line) starting from moderate num- 
bers of edge vacancies (see Figs. [8jb) for iV vac = 2 and 
gc) for N vac = 10). 

The features of P^'(S) presented above can be ratio- 
nalized with the ansatz for odd and even spacings 

P odd (a;S) = aP GOE (aS), (30) 
P cvcn (/3, k; S) = /3Pgoe-gue(k; /3S), (31) 

with P G oe(S), Pgoe-gue(A;S') given by Eqs. fl!5|23p 



and a constrain /? = a/ (2a — 1) guaranteeing, that the 
resulting distribution |[P dd + feven] satisfies (S) = 1 
for 1 a oo and K ^ oo. In Appendix |TJJ wc 
propose real random-matrix model with a single param- 
eter A controlling the transition from GUE with twofold 
level degeneracy (A = 0) to GOE without the degeneracy 
(A = oo) and utilize the ansatz ( 30|31 1. The empirical 
relations a = 5(A) (D4) and k = k(A) (D6 1 allow us to 



consider an approximating formula for nearest-neighbor 
spacing distributions 



p (l)» g s _ Podd(5(A);^) + P cvcn Q3(A),7;(A);S) 



For second-neighbor spacing distributions, we take 



(32) 



P™{\; S) = 2 dS'P odd (a(\); 2S - S') 
Jo 

i(/?(A), k(A); S'). (33) 



x P„, 



The empirical distributions P^(A; S) and P^%(X; S) 
with the parameter A = Afit (best-fitted for each value 
of A^vac) are shown in Fig. [8] with red solid lines. The 

(1 2) — 

asymptotic forms of (\;S) for A = and A = oo 

are depicted with blue dashed and blue dotted lines (re- 
spectively). In the first limit (A = 0) we have 



Pi%(0;S) = P GXJE (S), 



1 



S(S) 



i 



p 



GUE 



(S/2), 



(34) 
(35) 



restoring spectral properties of GUE with the exact 
twofold degeneracy of each level. The actual spacingdis- 
tributions for Ar vac = (see Fig.^a) for P«(S) and||d) 
for p( 2 )(S); black solid lines) show small deviations from 

(12) (12) 

(0; 5) and can be approximated by P±^ J (X^ t ;S) 
with Afit = 0.051, providing an estimation of the inter- 
valley scattering rate [75] in a triangular nanoflake with 
perfect zigzag edges. In the opposite limit (A = oo) [73] 



Pi%(^;S)cP GOE (S), (36) 

r2S 

p| 2 ^(oo; S) ~ 2 / dS'P GOE (2S - S / )P GO e(S / ) 
Jo 



7rexp (— ttS 2 



x exp 



irS 2 
2 



S 



erf 



irS 2 



V2 



S 



(37) 



The spacing distributions p( 1,2 )(S) for N vac = 10 [see 
Figs. Re) andliFf)] are close to P~ 2) (oo; S), but they 

— — . . (12) 

still fit to the approximating distributions (Afit; S) 

with Afit = 0.69 noticeably better. 

For the intermediate values of A^ vac , we generally ob- 
serve some systematic deviations of P^'(S) from the 
approximating distribution pl^Afit; 5) for S < 1 and 
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FIG. 9: (a) Least-squares fitted parameters Aat of em- 
pirical distributions P~L(\zt,;S) (|32[) approximating actual 

nearest-neighbor spacing distributions P^'(S) for triangular 
nanoflakes with zigzag edges and finite number of edge vacan- 
cies iV vac . Flake areas are A — (60 a) 2 (open symbols) and 
(120 a) 2 (close symbols), corresponding to the total number of 
terminal sites N c dgc = 270 and 540 (respectively). Solid line 
depicts the best-fitted power-law relation (381. (b) Spectral 
rigidity Ag(L) obtained numerically for A = (120 a) 2 with 
A^vac = (circles), Ar vac = 1 (diamonds), and iV vac = 10 (tri- 
angles). Lines correspond to theoretical expectations ( |28[ ) for 
GOE (red dashed) and for GUE with the twofold (valley) de- 
generacy of each level (black dotted). The disorder strength 
is fixed at Kq ~ 0.04 [7TJ for all systems. 
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FIG. 10: Level-spacing distributions P^\S) (a)-(c) and 
P (2) (S) (d)-(f) for the same system as in Fig. [sj (with N V3JC = 
0) placed in the uniform magnetic field (with the total flux 
<E> specified at each panel). Numerical results are shown with 
black solid lines. Other lines are for GUE with twofold de- 
generacy of each level, see Eqs. ( 34|35 [[blue dashed] ; and for 
two independent GUEs, see Eqs. ( |39|4o| [blue dotted]. 



good agreement for larger S. However, the comparison 
of PW(iS) shown in Fig. pi b) with the corresponding dis- 
tribution for smaller system [with the area A ~ (60 a) 2 ] 

suggests that PW (S) gradually converges to pl^Agt; S) 
with the system size also for S < 1. No significant devi- 
ations of P( 2 )(5) from P±%(\&u S) are observed. 

In Fig. |9ja) we plot the values of Agt for two different 
flake areas A ~ (60 a) 2 and (120 a) 2 (open and close sym- 

1/2 

bols, respectively) as functions of the variable N vac N e ^ ge , 
where A^dge denotes the total number of terminal sites. 
Such a scaling allows us to find a single power-law rela- 
tion for systems of different sizes, namely 



Am ~ 0.102(4) x 



N va . c N 



1/2 



edge 



0.34(1) 



(38) 



with the numerical values of parameters obtained via 
least-squares fitting (the standard deviation of a last digit 
are specified by numbers in parenthesis). The approxi- 
mating relation given by Eq. ( 38 1 is also depicted in Fig. 
[9]ja) [red solid line]. 

The evolution of more distant level spacings with in- 
creasing N vac is illustrates in Fig.^b), where we plot the 
spectral rigidity A 3 (L) for N vac 10. For N vac = 0, 
A 3 (L) obtained numerically closely follows the theoreti- 
cal expectation for GUE with the twofold valley degen- 
eracy of each level A 3 GUE) (L) = 4A 3 GUE) (L/2) Q (see 
circles and black dotted line, respectively). When in- 
creasing A/vac, the datapoints (diamonds or triangles for 
A/vac = 1 or 10) gradually approaches A 3 GOE ' ) (L) (red 



dashed line) for all values of L < 10. (With some devia- 
tions for larger L due to a finite system size.) 

Our demonstration of the nonstandard transition 
GUE-GOE in graphene nanoflakes, driven by varying the 
intervalley scattering rate at zero magnetic field, is now 
complete. Most remarkably, basic spectral characteris- 
tics start to reproduce those obtained in Refs. [T5H2H] for 
irregular edges, after removing just a few percent of ter- 
minal atoms from the system with perfect zigzag edges. 



B. Influence of external magnetic fields 

At finite magnetic fields, spectral statistics of graphene 
nanosystems with negligibly-weak intervalley scattering 
(a situation occurring for triangular nanoflakes with 
zigzag edges and N vac = 0) also require some atten- 
tion. In Fig. [lO] we plot level-spacing distributions for 
the flake area A ~ (120 a) 2 , and different magnetic fields 
B (quantified by the total flux piercing the system area 
$ = AB). The actual spacing distributions P^'(S) and 
P( 2 \S) (black solid lines), obtained numerically for the 
remaining system parameters same as in Fig. [HJ show 
crossovers from the theoretical predictions for GUE with 
twofo ld vall ey degeneracy (blue dashed lines) given by 
Eqs. ( 34|35 ) to the predictions for a level sequence follow- 
ing from two statistically-independent GUEs (blue dot- 
ted lines), approaching the latter for $ > $ - For such a 
combined sequence, we have the distribution of a Berry- 
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Robnik type [H] 



R 



(i) 



2xGUE 



(S) 



dS 2 



GUE 



(S/2)} 



(39) 



2 f 2S ! 

^2xGUe(^) = ^ / ^ -P2xGUe(25' — 'S' )^2xGUe(^' )> 

Jo 

(40) 

where ^gue(S') = e~ is '' V" - S + Seii{2S/^) [75] is 
the probability, that interval S contains no energy level 
of a simple sequence following GUE. 

The evolution of level-spacing distributions, presented 
in Fig. [lO] illustrates the fact, that for finite magnetic 
fields the effective Hamiltonian "HDirac ([l) does not com- 
mute with the symmetry T v Q and thus the valley-blocks 
are not degenerate. For finite systems, approximate val- 
ley degeneracy remains for $ <C $0 (as the valley energy 
splitting is much smaller than the average level spacing) . 
For <!> > <£>o> dynamical phases gained by carriers at K 
and K' valleys passing a typical closed trajectory start 
to differ significantly [76], and the sequences of energy 
levels corresponding to different valleys may be regarded 
as statistically- independent. 

It is worth to stress here, that the effect which we de- 
scribe may appear for real magnetic fields only. In con- 
trast, strain-induced pseudo-magnetic fields are exactly 
opposite at K and K' valleys, so they do not lift the val- 
ley degeneracy [131 [77] . We have found numerically, using 
the strained geometry of Ref . [78] , that level-spacing dis- 
tributions of chaotic nanosystems (each having the main 
symmetry axis bent into a piece of arc of the radius R) 
are unaffected even for extreme strains corresponding to 
R/H = 1. The details of the calculations will be pre- 
sented elsewhere. 



VI. SPECTRAL STATISTICS IN THE 
LOCALIZATION RANGE 

So far, the issue of the wavefunction localization in 
chaotic nanosystems in graphene has been addressed nu- 
merically in the literature by employing models of dis- 
order abruptly varying on the scale of atomic separation 
[79] . In this Section, we complement the existing studies 
with spectral statistics following from smooth impurity 
potential given by Eqs. (5|6|. The scope of the paper 
is thus extended on the examples of highly-symmetric 
graphene nanoflakes with strong potential disorder. 

It is know that in a generic quantum chaotic system 
eigenfunctions may not be uniformly distributed over a 
classically allowed phase space but localized (the dynam- 
ical localization effect) , which is associated with the frac- 
tional power-law level repulsion [80 . In the presence of 
TRS, one can expect the crossover from GOE in the case 
when extended chaotic states dominate the spectrum to 
the Poisson distribution in the strong localization limit 
|81j . The level spacing distribution for the system un- 
dergoing such a transition can be rationalized with the 



(b) K =0.99^ 
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FIG. 11: Nearest-neighbor spacing distribution P W (S) for 
triangular nanoflakes with zigzag edges, the area A = (120 a) 2 
[49] , and large values of the disorder strength Kq (specified 
at each panel). Numerical results are shown with black solid 
lines at all panels. The other lines at panel (a) and (b) are 
the same as in Figs.[8|a)-(c). Red solid lines at panels (c)-(f) 
correspond to the Brody distribution PBrody(/3; S) (411 with 
least-square fitted /3 = /3fl t ; blue lines are for the Poisson 
(dashed) and GOE (dotted) distributions. 
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FIG. 12: Least-squares fitted parameters /3fl t of Brody distri- 



butions PBrody(/3; S) (41 1 as functions of disorder strength for 
selected nanosystems of Table[T]in the localization range. Dat- 
apoints correspond to triangles with zigzag (A iVtot = 8278, 
TVtot = 32758) or armchair (T N tot = 32760) edges 
and smooth impurity potential (£ = ^/ia). The results 
for hexagons with zigzag edges and Anderson-type disorder 
(^ = 0) are also shown (Q Mot = 10584, • iV tot = 42336). 



well-known Brody distribution 

P Br ody(/3; S) = dSP exp(-C 2 S' 3+1 ), (41) 

with 

0+1 



C 1 = (0 + 1)C 2 , C 2 



13 + 2 



(42) 
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and T(x) being the Gamma function. The limiting distri- 
butions Pgoe(S) (15) or Pp ;(5) = exp(— S) are restored 
for /3 = 1 or /3 = (respectively). In turn, when ana- 
lyzing the spectral statistics of strongly disordered and 
closed nanosystem in graphene that preserves TRS, one 
can fit the empirical distribution Perody ■(/?; S) (41 ) to the 
actual nearest-neighbor spacing distribution P^>(S) m 
order to quantify the deviations from Pgoe(S), which 
indicate the localization. 

At sufficiently strong disorder all the systems studied 
in the paper show intervalley scattering which restores 
TRS (at zero magnetic field). In particular, the spac- 
ing distributions P^(S) obtained numerically for tri- 
angular nanoflakes with zigzag edges (N vac = 0), the 
area A = (120a) 2 , and different values of the disorder 
strength K [35] are depicted with black solid lines in 
Fig. [ill When increasing Kq, the distribution P^(S) 
first show a crossover from the theoretical prediction for 
GUE with twofold valley degeneracy [blue dashed lines 
in Figs.[llja) andJTlJb)] given by Eq. (34) to the Wigner 
surmise for GOE [blue dotted lines] , approaching the lat- 
ter at A'o — 1. We notice here, that for the system pa- 
rameters given in Ref. 49 Kq ~ 1 corresponds to the 
disorder amplitude S/t ~ 0.5, i.e., the orthogonal sym- 
metry manifests itself in P^(S) if \E\ < 5 for all en- 
ergy levels from the range 0.1 ^ \E\/t ^ 0.5, which are 
taken into account. Interestingly, in the crossover range 
(Kq < 1) PW(5) can still be rationalized with the empir- 
ical distribution P^l(X-S) (32) with least-square fitted 
A = Afit [red solid lines in Figs. 11 [a) andJTlJb)], sim- 
ilarly as in the case of transition GUE-GOE observed 
when increasing the number of edge vacancies N vac (see 
Sec.[v]). For K > 1, the distribution P^(S) can be ap- 
proximated by the Brody ditribution i-Brody (/?; S) (41) 
with least-squares fitted j3 — (3^ [red solid lines in Figs 
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c)-[TTTf )] and gradually approaches the Poisson distri- 
butionfblue dashed lines] signalling that the chaotic but 
localized eigenstates start to govern the spectrum. 
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The evolution of P^>(S) illustrated in Figs. 
|ll[ f) is qualitatively reproduced for all nanosystems con- 
sidered (see Table [i]) in their localization ranges. The 
value Kq, at which the crossover from GOE to the Pois- 
son distribution occurs, is however related to the sys- 
tem size and microscopic details of the disorder. For 
a quantitative description we plot, in Fig. [12] the val- 
ues of /3fit in the Brody distribution (41) best- fitted to 
the actual distributions P^'(S) obtained numerically for 
different nanosystems, as functions of the dimensionless 

1/2 

quantity iV tot Kq. After such a scaling of the indepen- 
dent variable, the available datasets group around just 
two distinct curves, one for the Anderson-type disorder 
characterized by £ = (circles in Fig. 12), and the other 



for smooth impurity potential with £ = V3 a (remaining 

1/2 

symbols). We attribute it to the fact, that iV to ' t K is 



proportional to the effective system size V Aj (1(E)) de- 
termined by the free path 1(E) in Born approximation 



1(E) 



2g v hvp 



1, if £ < a, 

2, if £ > a. 



Averaging 1(E) over the range E m [ n ^ |_Z?| -E-max; 



the weights Pbuik(^) given by Eq. (13), one gets 



2 • &l^g v t 



A 



E u 



■EL 



(1(E)) 



AAg v 



A 



(l(E)Y 



(43) 



with 



(44) 



where the last approximate equality holds true for 



E„ 



0.1 1 and E n 



0.5 1 used in our numerical 



simulations. 

It is also visible from Fig. [12] that the crossover to the 
localization range takes place for in case of £ = y/3 a for 

1/2 

iV tot more than two orders of magnitude larger than 
in the case of £ = 0. This suggests, that the smooth 
character of the disorder potential may be crucial for ex- 
perimental observation of signatures of quantum chaos in 
closed graphene nanosystems. 



VII. CONCLUSIONS 

We have investigated the symmetry classes of selected 
closed nanosystems in graphene (equilateral triangles 
and hexagons with three types of boundaries: armchair, 
zigzag, and Klein) and studied the effect of weak poten- 
tial disorder. Predictions of the Dirac equation for low- 
energy excitations were confronted with numerical results 
for the tight-binding model on a honeycomb lattice. New 
findings are visualized in Fig. [T] 

In the absence of disorder, available analytic solutions 
for continuous Dirac billiards show the level clustering 
due to number-theoretic degeneracies. Such degeneracies 
are lifted up in the tight-binding model due to nonlinear 
terms in the dispersion relation, leading to the Poissonian 
distribution of energy levels (characteristic for a generic 
intcgrable system). 

For weak disorders, the transition to quantum chaos is 
observed. In such a limit, spectral statistics follow these 
characterizing Gaussian ensembles of random matrices. 
In principle, all of the considered tight-binding Hamilto- 
nians are time-reversal invariant and expected to show 
the orthogonal symmetry class in the absence of mag- 
netic field. To the contrary, the Dirac Hamiltonian for 
graphene has a block structure related to the presence 
of two valleys making the true TRS irrelevant in the 
absence of intervalley scattering. Instead, special TRS 
(symplectic symmetry) applies and is broken by the dis- 
order, leading to the unitary class (accompanied by the 
twofold valley degeneracy of each level). In effect, the 
type of boundaries plays a decisive role for the symmetry 
properties. 

Earlier studies of closed nanosystems in graphene [18l - 
[2171 12"4"] reported the orthogonal symmetry class associ- 
ated with the valley mixing. We have found that the 
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unitary symmetry class can also be observed in spec- 
tral statistics of such systems, providing almost all termi- 
nal atoms belong to one sublattice. This is satisfied for 
equilateral triangles with zigzag or Klein boundaries, for 
which the spectral statistics obtained numerically show 
the following features: When increasing the disorder 
strength, transition from the Poisson to GUE distribu- 
tion (both showing an approximate twofold degeneracy 
of each level) occurs. For a fixed disorder strength in the 
chaotic range and increasing the number of edge vacan- 
cies we have observed the transition to GOE distribution 
(accompanied by the gradual level splitting). Moreover, 
for the same disorder strength and in the absence of edge 
vacancies, we have demonstrated the evolution from GUE 
distribution with twofold degeneracies to the distribution 
characterizing two independent GUEs at weak magnetic 
fields. These findings complement the very recent results 
[23] for transport characteristics of open nanosystems in 
graphene. 

The remaining nanosystems studied in the paper are in 
the orthogonal symmetry class. For all cases, the tran- 
sition to quantum chaos is rationalized using additive 
random-matrix models. The functional relation between 
the best-fitted model parameter Afi t and the disorder 
strength has a form of a power law Agt oc (N to tKo) a , 
with the symmetry-dependent exponent a taking differ- 
ent values for systems undergoing the transitions Poisson- 
GOE and Poisson-GUE. Additionally, the model involv- 
ing 4x4 real random matrices is proposed and elabo- 
rated to parametrize the nonstandard GUE-GOE tran- 
sition identified for triangular flakes with edge vacan- 
cies at zero magnetic field. For strong disorders (i.e., 
in the localization range) additive random-matrix mod- 
els no longer apply. Instead, the fractional level repulsion 
and the evolution, with the increasing disorder strength, 
towards the Poissonian distribution of energy levels are 
observed, indicating the spacial localization of quantum 
states. 

We hope recent progress in resolving closely-lying en- 
ergy levels in graphene quantum dots using the three- 
terminal Coulomb-blockade setup |83| will make it pos- 
sible to test experimentally our results. 
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Appendix A: Boundary conditions to Dirac equation 
and symmetry classes of graphene nanoflakes 

Boundary conditions for Dirac fermions in graphene 
are usually discussed in the so-called valley-isotropic rep- 
resentation [HE! • I n this Appendix, we recall the standard 
expressions (3] [30] for infinite mass, armchair, zigzag and 
Klein boundaries, and rewrite them in the notation of 
Eq. ([I]) to illustrate how particular boundaries may de- 
termine the system symmetry class. 

The valley-isotropic representation of the Hamiltonian 
([I]) can be obtained using the following unitary transfor- 
mation [m EI] 

"H D irac = W t- HDirac U = V F [(p + eA) ■ <t] <g> T 

+ M(r)<r z (g) t z + U(r)a ® r , (Al) 

where 



U 











o\ 





1 
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\o 





-1 





2 CT o ® (to + T z ) + -o-y (t 



r z ), (A2) 



and the remaining sym bols are the same as in Eq. ([I]). 
The Hamiltonian (Al) now acts on spinors ip = U'il> — 
[ipA^B, —tpB^AV ■ For model situations considered in 
the literature [301 [47] the mass term M(r) = and the 
Hamiltonian (Al) contains only terms proportional to 
To, so it consists of two identical blocks (one for each 
valley) justifying the notion of 'valley-isotropic represen- 
tation'. As we show in Section [lI B[ the potential disorder 
in graphene leads to M(r) ^ 0, so the term proportional 
to t z appears in the Hamiltonian. However, the rep- 



resentation (Al A2) still remains useful for defining the 
boundary conditions to Dirac equation. Also, the current 
operator 



vper <g> tq 



(A3) 



is proportional to a ® To and thus has identical form for 
both valleys, regardless M(r) = or M(r) ^ [So], 

Most common boundary conditions for graphene 
nanosystems may be expressed in a compact form as [85] 



ijj = Mip, M = (rj ■ er) <8 (v ■ r) , 



(A4) 



where rj and v are three-dimensional unit vectors. The 
vector r\ is constrained by r\ ■ r\ B to guarantee that no 
current leaks out of the boundary, defined by the nor- 
mal r\ B (pointing outward). In fact, Eq. (A4) represents 



the general boundary condition, providing we ignore non- 
collincar local magnetization (which may appear on the 
edges of a graphene nanoflake [57] [EE]), so one can as- 
sume that the boundary condition itself does not break 
time-reversal symmetry [3]. 

The examples of boundaries that can be defined via 
vectors r\ and v in Eq. (A4) are: 
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• Confinement by an infinite mass corresponds to 
rj = ± e z x r\ B and v = e z , where the upper (lower) 
sign is valid for the mass going to +00 (—00) out- 
side the system area. 

• An armchair edge requires the wavefunction 
([2]) is vanishing on both sublattices, namely: 
ipx exp (iK ■ r) + ip' x exp (— iK ■ r) = for X = 
A, B. This corresponds to r\ = ± e z x r] B and 
e z • v = 0, where the upper (lower) sign is valid 
when the order of the atoms within each dimer is 
A—B (B — A) along the direction of e z x rj B . 

• A zigzag (or Klein) edge requires ipA — ^P'a — or 
ipB — ip's = 0, depending on whether the row of 
missing atoms at the edge is on the A or B sublat- 
tice, what corresponds to 77 = ± e z and v = e z . 



The boundary condition (A4) can be written in the 
notation of Eq. ([!]) as 



(A5) 




<fA 
B 



where 



FIG. 13: Left — two zigzag edges of different types forming a 
120° corner, right — an armchair edge (tick solid lines). These 
two boundaries become equivalent if an electron (or hole) in- 
coming from the sample area (shaded) has a Fermi wavelength 
\f large compares to the lattice spacing a. The atoms in 
terminal and first missing rows are indicated with open or 
closed dots depending if they belong to A or B sublattice. 
(Each dashed line marks a bond between the terminal and 
the nearest missing atom.) r) B , r)' B , and rfg are unit vectors 
normal to the boundaries. The coordinate system and the 
first Brillouin zone (inset) are also shown. 



M = {r) x a x + r) y (jy) <X> v z t q + r\ y a y ® v z t z 
+ (Vx&z - n z a x ) <E> (y x r x + v y T y ) 

+ %cr <g> {v x T y - v y T x ) . (A6) 

For the most common boundary conditions listed above, 
M is given by 



M = 



1] X <T X ®T + T) y Oy <g> T z 

T] X G Z ® (U X T X + UyTy) 

+ %CT <g> (V X Ty-VyT X 

±O z ®Tq 



(infinite mass) 

(armchair) 
(zigzag/Klein) 



(A7) 



where r\ — (r\ x -,r\y) — ±e z x r] B for the first two cases. 
For armchair edge, v = {y x ,Vy) is a unit vector in the 
x — y plane. 



It is clear from Eq. (A7) that the boundary condition 



( A5 ) couples the valley degree of freedom for the case of 



armchair edges, leading to the orthogonal symmetry class 
of a chaotic nanosystem at zero magnetic field. For the 
remaining cases, no obvious intervalley scattering origi- 
nates from the edges, so one may expect that the unitary 
symmetry class appears. The more carefully discussion 
is necessary, however, for systems containing two zigzag 
(or Klein) edges of different types. 

In particular, we focus here on 120° corners formed 
by a zigzag edge terminated on A sublattice attached 
to a zigzag edge terminated on B sublattice (such as 
presented in Fig. 13), which appear in hexagons with 



zigzag edges. An effective wavefunction for low-energy 
excitations ^ near the lower arm of such a corner is 
determined by the boundary condition that may be writ- 
ten as ipA exp (iK • r) + ip' A exp (— iK • r) = 0, where r 



denotes the position at the first row of missing atoms. 
This condition cannot be satisfied simply by setting 



0, because the wavefunction also needs to 



satisfy an analogous condition for the upper arm, namely: 
tpB exp (iK ■ r') + ip' B exp (— iK ■ r') = 0. In the coordi- 
nate system of Fig. [l3j we immediately get that these 
two conditions lead to 



V 



e z x rf§ 



V C B S 



e z v 

Vb + Vb 
\Vb+V'bW 



0, 



(A8) 



This is an effective armchair boundary condition. 

We can expect now, that strong intervalley scattering 
originating from 120° corners places hexagons with zigzag 
edges in the orthogonal symmetry class. Similar argu- 
ments apply to hexagons with Klein edges. For this rea- 
son, the only closed systems for which the unitary sym- 
metry class may still manifests itself in spectral statistics 
are those bounded entirely with zigzag (or Klein) edges 
with terminal atoms belonging to one sublattice, such as 
equilateral triangles considered in Sections |IIIffV] 



Appendix B: Low-energy level structure of 
triangular graphene flakes 



In Tables 



i|ffl|and|rv|^ 

iguiar nanotle 



I we list energy levels -E°% a + (10) 
of two triangular nanonakes considered in the paper (see 



Table |T|, by taking n < 5 and m-s satisfying ( 11 ) or ( 12 ) 
for armchair or zigzag edges (respectively). Correspond- 
ing energies obtained from the exact numerical diagonal- 
ization of tight-binding Hamiltonians ^ in the absence 
of disorder (Kg = 0) and for infinitesimally-weak bulk 
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TABLE III: Lowest-lying electronic levels (JJ°™ + ) obtained 
from Eq. ( 10 1 for triangular graphene flake with armchair 



TABLE IV: Same as Table [fflj but for triangular flake with 
zigzag edges and iVtot = 32758. 



edges containing Aftot = 32760 atoms and corresponding en- 
ergies (E) obtained from numerical diagonalization of Htba 
{5} for zero and for weak disorder (specified by Kq). The 
staggered potential Mv{ri) = for all lattice sites. 
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pDirac l± 


E/t, 


Ki 
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-fiy i n -4 
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0.02004229 


0.01994629 





01992564 






0.01994629 





01992624 
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0.03474520 
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0.03989059 





03960217 




u.uoouzuyi 


U.UJiJ / oou 


n 








0.05315689 





05271721 


(1,3) 


0.05302691 


0.05237556 





05225943 






0.05315689 





05312196 




U.UDUlZDoO 


u.uoyoouyo 


u 








0.05983093 





05990582 






n nR898^i r 

U.UUOZOOIU 


u 


flfi7Q7fi91 






0.06987442 





06954428 


(3,4) 


0.07226350 


0.07132749 





07112046 






0.07247406 





07228640 




U.Ui zzuoou 




n 

u 


071 ^RQI Q 






0.07247406 





07272521 


(4,4) 


0.08016915 


0.07976532 





07971228 






0.07976532 





07971380 


(2,5) 


0.08736231 


0.08572437 





08523396 






0.08809530 





08767157 


(3,5) 


0.08736231 


0.08572437 





08566093 






0.08809530 





08811008 


(4,5) 


0.09184530 


0.09062032 





09035237 






0.09212365 





09184820 


(1,5) 


0.09184530 


0.09062032 





09059671 






0.09212365 





09209424 


(5,5) 


0.10021144 


0.09969177 





09933460 






0.09969177 





09936139 



(m, n) 


pDirac /j. 
- C/ m,n,4-/ 1 


E/t, 
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E/t, 
A"o = 6xlO- 4 


(1,2) 


0.03471532 


0.03452128 
0.03452128 


0.03437910 
0.03437919 


(1,3) 


0.05302853 


0.05247781 
0.05247781 
0.05297828 
0.05297828 


0.05208748 
0.05208754 
0.05301960 
0.05301991 


(2,4) 


0.06943064 


0.06903228 
0.06903228 


0.06873980 
0.06874017 


(1,4) 


0.07226570 


0.07120437 
0.07120437 
0.07248965 
0.07248965 


0.07120256 
0.07120295 
0.07228947 
0.07228985 


(2,5) 


0.08736497 


0.08642740 
0.08642740 
0.08727781 
0.08727781 


0.08616840 
0.08616900 
0.08701618 
0.08701620 


(1,5) 


0.09184810 


0.09012997 
0.09012997 
0.09247007 
0.09247007 


0.09000825 
0.09000841 
0.09214693 
0.09214700 



disorder (Kq ~ 10~ 3 ; for the disorder details see Ref. 
|52j ) are also provided. Unlike for similar Schrodinger 
systems [Ml 84] , the spinor structure of the wavefunction 
([2]) cause that geometric symmetries of graphene flakes 
do not lead directly to level degeneracies. Instead, degen- 
eracies associated with special time-reversal symmetries 
%i and T v ^ may appear. 

We see from Table IIIIl that for armchair boundaries 



each electronic level of the Dirac cavity E^ I ^ C , is followed 
by two levels of the lattice system, corresponding to the 
Kramer's degeneracy associated with %i Q. For the lat- 
tice system, this degeneracy is usually only approximate 
even at Kq = 0, as a nonlinear term appearing in the 
effective Hamiltonian derived from tight-binding model 
does not commute with 7l\ |35j . Analogous symmetry T v 
Q does not apply for armchair edges, however, the prop- 
erty Em™ c = En-m,n lea ds to an additional twofold de- 
generacy of each level, providing that m ^ n and 2m n. 
As a result, for Kq = majority of energy levels occurs in 
almost-degenerated quadruplets. The degeneracy is im- 
mediately lifted in the presence of disorder. Typically, as 
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small disorder strength as Kq ~ 10 3 leads to the level 
splitting SE > 10~ 4 i for \E\ < 0.1 1. 

For zigzag edges and for My = (see Table IV), the 



fourfold (approximate for the finite lattice system) de- 
generacy appears for almost every level, due to the sym- 
metries 7^i and T v Q- For even n, the degeneracy is only 
twofold when 2m = n. We attribute this to the presence 
of edges states in the system with zigzag edges, which 
have missing valley degeneracies of bulk states. Such 
levels, however, do not contribute to spectral statistics 
of large systems. Unlike for armchair boundaries, the 
twofold valley degeneracy associated with the symmetry 
T v Q is present for almost every level of the triangle 
with zigzag edges, and appears to be very robust against 
the disorder. For K ~ 10~ 3 , corresponding splittings 

< 



from Table 



-3 

IV] are SE < 10- 7 'i for \E\ < 0.1 1. (Notice 
that the degeneracy associated with T s \ is lifted in the 
presence of disorder for either zigzag or armchair bound- 
ary conditions.) These splittings are unaffected when the 
staggered potential with |My(rj)| = 0.7 1 is put on the 
outermost edge atoms. The identical structure of energy 
levels as for zigzag edges was observed for the case of a 
triangle with Klein edges. 

These are the reasons, why large and weakly- 
disordered triangular nanoflake in graphene shows the 
twofold, approximate level degeneracy only if it has per- 
fect zigzag (or Klein) edges. In other cases, no degenera- 
cies appear in the presence of disorder. 



Appendix C: Level clustering in triangular graphene 
flakes 



Electronic energies of triangular Dirac cavities ( 10 ) 
may be written, in the dimensionless units, as 



= 2^fr 



W- — mn = 



y/(m + n) 2 + 3(to - n) 2 = ^Jl 2 + 3k 2 = e kl , (CI) 



where A = irt/\/3N to t, k = n — to, and I = n + to 
[55] , Without the loss of generality, we have limited the 
discussion to energy levels in the conduction band. Let 
tki ^ e^'i' are the nearest neighbors in the sequence. For 
high energies, i.e., ekli^k'V 3> A, we have 

e 2 kl - <? vv = l 2 + 3k 2 {l'f - 3{k'f - 

2ejw (ekl - ejfe'i') = ^^-p(eki) (e k i - £k>i>) , (C2) 

7T 

where the dimensionless density of states p(e) is defined 
via Pbuik(E) = A^ftE/A), with p huik {E) given by Eq. 
(fl3b- By definition, p(e k i) (e k i - £k'l') = S. Denoting the 



integer in the first line of Eq. ( C2 ) by 

q = I 2 + 3k 2 - (I') 2 - 3(k') 2 



(C3) 



we immediately obtain the quantization rule for S given 
by Eq. (17) in the main text. 




FIG. 14: Average level degeneracy for a triangular 
nanoflake with armchair edges containing A to t atoms for en- 

and 



ergies J? m',T?,± o btained from Eq. (10 1 [blue solid line 

( C4|C5 l [black dotted line] all taken from the range 



( — 1/2, i/2j! The asymptotic form g n b,o (C6l is also shown 
[red dashed line]. Inset shows the deviation of the actual de- 
generacy from <7nb,o in a magnified vertical scale. 



Similarly as for the analogous Schrodinger systems 
55, 56, 84 , energy levels of equilateral triangles contain- 
ing Dirac fermions show number-theoretic degeneracies 
not connected to the geometric symmetry. In particu- 
lar, the result by Pinsky [56] who showed that the av- 
erage level multiplicity is divergent when expanding the 
energy interval from which the levels are taken into ac- 
count applies directly to energy levels given by Eq. (10). 



The divergence (or level clustering) also appears when 
fixing the energy range, i.e., \E^"^ ± \ E max , and in- 
creasing N tot . For these reasons, equilateral triangles 
with Dirac fermions cannot be regarded as generic inte- 
grable systems, unless nonlinear terms in the dispersion 
relation (originating from the tight-binding Hamiltonian 
of graphene) lift up number-theoretic degeneracies. 

To further illustrate a possible role of number-theoretic 
degeneracies in graphene nanoflakes we consider the en- 
ergy levels recently found by Rozhkov and Nori for a 
tight-binding Hamiltonian of the equilateral triangle with 
armchair edges [4"5] 



pTBA 



+ 2 COS 



±t 



^3 + 2 cos 


27m 


.3(AT a + l). 





27TTO 

3(JV a + l) 



2 cos 



27r(n 



3(JV„+1) 



1/2 



(C4) 



where N a = 3H/(2a) (such that AT tot = 3N a {N a + 1), 
see Ref. [H]), fa = N a — to, and h = N a + n. For 1 ^ 
to n -C N a one gets 3(iV Q + 1) ~ y/3N tot and #™A ± ~ 
^ron± ' res t° rm g the energy l evels of a Dirac cavity ( 10 1 . 



Furthermore, expanding Eq. (C4| in series and keeping 
the terms up to the order of ~ m r n^ r ^ / {N a ) 3 , with the 
integer r ^ 3, we can write 



pTBA 

m,n,± 



(efc,i) - 



_tt_ I 2 (I - k) 
~j% N a + 1 : 



(C5) 
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where the notation of Eq. (CI) is used. The last term in 



Eq. ( C5 ) represents so-called trigonal warping of the dis- 



persion relation appearing in the vicinity of each Dirac 
points in graphene [35] . It is clear, that the quantity 



( pTBA 1 



(-^™«'.±) 2 i s n °t an integer multiplicity of 



A for arbitrary (m',n') ^ (m,n), and thus the quanti- 



zation rule for S (17) no longer applies 



Additionally, the trigonal warping appears to be the 
reason for which number-theoretic degeneracies are to- 
tally absent in spectra of tight-binding Hamiltonians 
for graphene nanoflakes. In Fig. [l4j we plot the aver- 
age degeneracy (<? n b) of energy levels from the interval 
(— t/2, t /2) obtained from Eqs. (10) [blue solid line], and 
( |C4|C5| ) [black dotted line] as a function of iVtot ^ 10 . 
(Notice that as we focus on number-theoretic degenera- 
cies, the twofold degeneracy = P™„ „.± of each 
energy level for which m^nis not taken into account.) 
The asymptotic form of g n b for large N to t in the presence 
of number-theoretic degeneracies |55j 



5nb,o oc vAog (TVtot /const) 



(C6) 



is divergent for AT to t — > oo and depicted in Fig. [M] [red 
dashed line]. For Dirac cavities, energies E^£ c ± (10) 
show the first number-theoretic degeneracy at -/Vtot = 
2610, for which g nh = 1| > 1. For larger AT tot , g nh ~ g nh 
with the accuracy better than 1% if N to t ^ 10 6 . In con- 
trast, neither exact energies of tig ht-b inding Hamiltonian 
for triangular nanoflakes E^^ ± ( C4 ) nor these given by 
the approximating Eq. ( C5 ) show number-theoretic de- 



generacies (i.e., g nh — 1 for any N tot ). 



Appendix D: Transition GUE-GOE for 4 x 4 real 
symmetric matrices 



a = 6.405(3), K = 0.93(2) 




FIG. 15: Level-spacing distributions averaged over 10 
randomly-generated matrices H(X) [datapoints] . The scaling 
parameter A is varied b etwe en the panels. The least-square 



fitted functions P^ K (S) (D3l are also shown [solid lines 



which has a twofold eigenvalue degeneracy and belongs 
to GUE. The matrix V in Eq. ( 19 ) is now chosen as a 
2N x 27V member of GOE. 

Varying the scaling parameter A, we transform the 
symmetry class of random matrix H(X) from the uni- 
tary (A = 0) to the orthogonal (A = oo), simultaneously 
splitting the twofold eigenvalue degeneracy. The nearest- 
neighbor spacing distributions are approximated by 



-Podd(a; S) + P C vcn(/3, k; S) 



(D3) 



In this Appendix, we analyze numerically a simple 
additive model of random matrices capable of describ- 
ing the transition GUE-GOE associated with the split- 
ting of twofold valey degeneracy in graphene. A single- 
parameter formula approximating the nearest-neighbor 
spacing distributions PW(5) is proposed. 

The analysis starts from the model 2N x 2N real sym- 
metric matrix H(A) of the form (19), where H° has a 
block-structure 



H° = 




(Dl) 



with A = A T and B = —B T . The elements of each block 
are independently generated according to a Gaussian dis- 
tribution with zero mean and the variances Var(Ajj) = 
(l+%)/2ATandVar(-Bi?) = (l-£y)/2iV. Therefore, H° 
can be unitary mapped onto complex Hermitian matrix 



H° = 



A + iB 
A + iB 



(D2) 



with P dd(o!; S) and P eV e n(/3, k; S) given by Eqs. ( 30 ) and 



(31), respectively. Eq. (D3) represents an observation, 



that the spacings distribution for the random matrix 
H(X) consists of two contributions (both of the equal 
weights for large N): first from odd spacings, separating 
the levels that are degenerate for A = 0, and the second 
from even spacings. We further suppose the distribution 
of odd spacings is well approximate by the Wigner sur- 
mise for GOE (15), whereas even spacings undergoes the 



transition GUE-GOE according to the Berry-Robnik for- 
mula ( 23 1 . The relation between parameters 1 ^ a < oo 



and ^ k < oo of the spacing distribution P Q , jK (5) and 
the scaling parameter A is to be determined. 

We test numerically the formula (D3) for N — 2 and 
A = 10~ 2 — 10 2 . For each value of A, an ensemble of 
10 7 — 10 8 pseudorandom matrices was generated, each 
4x4 matrix was diagonalized, and a histogram of spac- 
ings distribution pW(S') was obtained. To truncate the 
large- N limit, in which the contributions from odd and 
even spacings have equal weights, we took each first and 
third spacing with the weight 1/4, whereas second spac- 
ing weight was set to 1/2. Subsequently, the function 
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a = 1.118, K = 




FIG. 16: Deviation of the spacing distribution for 10 ma- 
trices H(X) with A = 10 from the Winn er surmise for GOE 
Jlj|. Blue solid line shows P£l(S) ( D3 1 fitted to the actual 



data [points]. Red dashed line corresponds to the asymptotic 
(N -> oo) distribution for GOE [90]. 



lows us to propose empirical functions a = a (A) and 
K = k(A). The first function reads 



a(A) = a 



V 









i 


i+ 











b/c 



(D4) 



with 



a = 1.118(1), 



Ai = 0.607(3), 6 = 0.978(2), 
c= 3.2(4), (D5) 



where standard deviations obtained from least-square fit- 
ting are specified in the parenthesis. (It is worth to stress, 
that the value of a parameter ag was found directly from 
spacing distributions corresponding to A > 10 for which 
k ~ 0.) The second function is given by 




1+X~ 



Best-fitted parameters of Pi]l (S) ( D3 1 [datapoints] 



FIG. 17: 

and the empirical functions a 

[solid lines]. Dashed lines depict the asymptotic forms 



a(A) (|D4| and n = k(A)JD6| 
D8T 



Po,k{S) was fitted to the numerical data within the least- 
squares method. Selected examples are presented in Fig. 
To"] Remarkably, the parameter a does not approach 1 
for large A. This is because Wigner surmise Pgoe(S) 
( 15 1 represents the spacing distribution exactly for 2x2 



GOE matrices only. For 4x4 matrices, formula (|D3j) with 
a = ao ~ 1.118 and n = appears to provide much bet- 
ter approximation of actual spacings distribution P(S) 
than the Wigner surmise (see Fig. 16). 

The collection of best-fitted distributions Pi]l(S) al- 




(D6) 



with 



7o = 0.2898(5), A c = 0.329(2), 6 = 2.6(1), 

A 2 = 0.178(2). (D7) 



The functions 5(A) and k = «(A) are plot in Fig. 17 
(solid lines) together with actual datapoints used for the 
fitting. The asymptotic expressions for small A 



a (Ai/A) 6 



1 + A- 
l + A; 
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(D8) 



are also shown in Fig. nj\ (dashed lines). Substituting 
a = 5(A) and n — k(A) to the formula (D3| we ob tain 
a single- parameter function pl x l(A; S) given by Eq. (32) 
in the main text. Thus, the construction of an empirical 
formula for the nearest-neighbor spacing distribution of 
H(X) is complete. 
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